//==============================================================================
//3456789_123456789_123456789_123456789_123456789_123456789_123456789_123456789_
//                               (why all the numbers? I prefer 80-column text)
//
// CPartSys.cpp         // Member function bodies for the CPartSys class.
//==============================================================================

#include "CPartSys.h"   // the *ONLY* #include statement in this file.
                        // (if you need to #include another file,
                        //  put the #include statement in CPartSys.h, not here!)

CPartSys::CPartSys(void)
//------------------------------------------------------------------------------
// Default constructor
{
    pS0 = pS0dot = NULL;        // make state-vector pointers benign.
    pS1 = pS1dot = NULL;
    pSM = pSMdot = NULL;
    pF0 = NULL;                 // and forcer-vector pointer benign,
    pC0 = NULL;                 // and wall-vector pointer benign.
    partCount = 0;              // and set the number of particles to zero.
    forcerCount = 0;
    wallCount = 0;
    runNotStep = false;             // default: start in single-stepping mode.
    timeStep = 0.01;                // default timestep
    solverType= NU_SOLV_EULER;      // default; use Euler or Forward solver.
}

CPartSys::~CPartSys(void)
//------------------------------------------------------------------------------
// Default destructor
{
    stateVecWipe(); // discard all the particle objects (CPart) we allocated for
                    // all of our state vectors pS0, pS0dot, pS1, etc.
   forcerVecWipe(); // discard all the force-making objects (CForcer) we alloc'd
                    // for all our forcer-vectors pF0, etc.
   wallVecWipe();   // discard all the constraint-making objects (CWall)
                    // that we alloc'd for all our limit-vectors pC0
    return;
}

void CPartSys::stateVecSizer(int nuSize)
//------------------------------------------------------------------------------
// If state vectors currently have any particles, DELETE them all; make empty
// state vectors.  Then for each state vector create a new dynamically alloc'd
// array that holds 'partCount' particles (CPart objects) and point to it.
//(remember, state vector is a ptr. to the first in an array of CPart objects.)
{
    // alloc a new array of CPart objects, one array for each state variable;
    if(nuSize >0)               // For any non-zero size;
    {
        if(pS0 != NULL) stateVecWipe(); // discard any existing state vectors.
        pS0     = new CPart [nuSize];   pS0dot  = new CPart [nuSize];
        pS1     = new CPart [nuSize];   pS1dot  = new CPart [nuSize];
        pSM     = new CPart [nuSize];   pSMdot  = new CPart [nuSize];

        if( pS0==NULL || pS0dot==NULL // if ANY allocation failed,
        ||  pS1==NULL || pS1dot==NULL
        ||  pSM==NULL || pSMdot==NULL)
        {
        cout << "\n!!!CPartSys::stateVecSizer() allocation error!!!\n" << endl;
        }
        else                            // if it succeded, record array size;
        {
            partCount = nuSize;
        }
    }
    else
    {
        cout << "!?! CPartSys::stateVecSizer(val) val<=0 ?!? call stateVecWipe() instead!" << endl;
        stateVecWipe();
    }
}

void CPartSys::stateVecWipe(void)
//------------------------------------------------------------------------------
// Discard all existing memory allocated for state vectors in this CPartSys obj.
{
    if(pS0!=NULL || pS0dot!=NULL ||
       pS1!=NULL || pS1dot!=NULL ||
       pSM!=NULL || pSMdot!=NULL)   // Do we have *ANY* CPart objects ANYWHERE?
    {   // Yes.  Then we should have the same number of CPart objects at each
        // and every pointer.  Delete them all, but complain if any are missing:
        if(pS0==NULL) // make sure this vector isn't empty, and
            cout << "CPartSys::stateVecWipe()--empty s0 vector ?!?" << endl;
        else delete(pS0);   // empty it.
        if(pS0dot==NULL)
            cout << " CPartSys::stateVecWipe()--empty s0dot vector ?!?" << endl;
        else delete(pS0dot);
        if(pS1==NULL)
            cout << " CPartSys:: stateVecWipe()--empty s1 vector ?!?" << endl;
        else delete(pS1);
        if(pS1dot==NULL)
            cout << " CPartSys:: stateVecWipe()--empty s1dot vector ?!?" <<endl;
        else delete(pS1dot);
        if(pSM==NULL)
            cout << " CPartSys:: stateVecWipe()--empty sM vector? !?!" << endl;
        else delete(pSM);
        if(pSMdot==NULL)
            cout << " CPartSys:: stateVecWipe()--empty sMdot vector ?!?" << endl;
        else delete(pSMdot);
    }
    else    // No... ?? Why you call me with an empty vector!?!?
    {
        cout << "\nCPartSys::stateVecWipe() called with all stateVec objects empty?!?\n" << endl;
    }
}

void CPartSys::stateVecZero(CPart *pS)
//------------------------------------------------------------------------------
// Set to zero all members of all particles in the state vector at pS.
// Useful for initializing state-variable derivatives; if all are zero, then
// by default we expect NOTHING to change over time, then set the exceptions.
{
int i,j;

    if (pS==NULL || partCount==0)
    {
        cout << "\n!!!ERROR!!! CPartSys::stateVecZero() called on zero-sized state vector!!?!\n" << endl;
        return;
    }
    for(j=0; j<partCount; j++)      // for every particle in the state vector
    {
        for(i=0; i< PART_MAXVAR; i++)  // for every variable in that particle,
        {
            pS[j].val[i] = 0.0;      // set it to zero.
        }
    }
}

void CPartSys::stateVecAplusBtimesC(CPart *pDest,
                                    CPart *pA, CPart *pB, double sc)
//------------------------------------------------------------------------------
// Compute with state vectors:  replace pDest state vector contents with
// A + sc*B;
// REQUIRES that the size of A and B are both nonzero and match exactly.
// Exactly what we need for Euler integration!  Far more elegant: use C++
// operator overloading to make your own versions of =, +, -, * for state vars.
{
int i,j;
    // error checks first:
    if(NULL==pDest || NULL==pA || NULL==pB || partCount <= 0)
    {
        cout << "\n!!!CPartSys::stateVecAplusBtimesC() FAILED: empty state vec!?!?\n" << endl;
        return;
    }
    for(j=0; j< partCount; j++) // For each particle in our state vectors,
    {
        pDest[j].partType = pA[j].partType;
                                        // copy 'type' from A (usually s0)
        for(i=0; i< PART_MAXVAR; i++)   // for each variable in that particle,
        {   // replace our value with A + sc*B
            pDest[j].val[i] = pA[j].val[i] + sc*(pB[j].val[i]);     // copy.
        }
    }
}

void CPartSys::stateVecAssign(CPart *pDest, CPart *pSrc)
//------------------------------------------------------------------------------
// copy contents of state vector at 'pSrc' to state vectir at 'pDest'.
{
int i,j;

    if(NULL==pDest || NULL==pSrc || partCount <= 0 )
    {
        cout << "\n!!!::CPartSys::stateAssign() FAILED: size mismatched !?!?\n" << endl;
        return;
    }

    for(j=0; j< partCount; j++)          // For each particle we have,
    {
        pDest[j].partType = pSrc[j].partType;   // copy pesky 'partType' flag
        for(i=0; i< PART_MAXVAR; i++)   // for each variable in that particle,
        {
            pDest[j].val[i] = pSrc[j].val[i];     // copy.
        }
    }
}

void CPartSys::stateVecSwap(CPart **ppA, CPart **ppB)
//------------------------------------------------------------------------------
// Swap the contents of state vectors pA, pB (which each are a pointer-to-CPart)
// TRICKY!! if you want to swap state vectors pS0 and pS1, you must call this
// fcn with the ADDRESS OF pS0 and pS1, like this:
// stateSwap(&pS0, &pS1);  Then ppA[0] holds pS0's value, ppB[0] holds pS1's.
{
CPart *pTmp;

    if(0==partCount || NULL==ppA[0] || NULL==ppB[0])
    {
        cout << "\nERROR!!! CPartSys::stateSwap() mismatched vectors!?!?\n" << endl;
        return;             // error exit
    }
    // Swap pointer values:
    pTmp = ppA[0];              // save A
    ppA[0] = ppB[0];            // A gets B
    ppB[0] = pTmp;              // B gets saved A.
}

void CPartSys::stateVecDraw(CPart *pS)
//------------------------------------------------------------------------------
// Draw on-screen all CPart objects found at the state-vector pointer pS.
// (number of CPart objects is CPartSys member 'partCount'
{
int i;
GLfloat P_diffuse[]   = {1.0, 0.5, 0.1, 1.0};	// gold,
GLfloat P_specular[]  = {1.0, 0.8, 0.4, 1.0};  // gold specular,
GLfloat P_shiny = 50.0;

	glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, P_diffuse);
	glMaterialfv(GL_FRONT, GL_SPECULAR, P_specular);
	glMaterialfv(GL_FRONT, GL_SHININESS, &P_shiny);

    // in whatever openGL coordinate system we're currently using,
    // in whatever glMatrixMode() settings we have (should be GL_MODELVIEW)
    // Draw each particle;
    for(i=0; i< partCount; i++)         // for each particle,
    {
        if(pS[i].partType <= 0) continue;         // skip dead particles
        glPushMatrix();                     // save current matrix, then
            glTranslated(pS[i].val[PART_XPOS],    // translate to position
                         pS[i].val[PART_YPOS],
                         pS[i].val[PART_ZPOS]);
            glutSolidSphere(0.04,7,5);     // Draw the particle as a sphere,
        glPopMatrix();                      // restore the saved matrix.
    }
}

void CPartSys::forcerVecSizer(int fsiz)
//------------------------------------------------------------------------------
// if forcer-vectors currently have any CForcer objects, DELETE them all; make
// empty forcer-vectors.  Then for each forcer-vector, create a new dynamically
// allocated array that holds 'forcerCount' forcers (CForcer objects), and
// point to it.
// (Remember, forcer-vector is aptr. to the first in an array of CForcer objs.).
{
    if(fsiz > 0)                        // for any non-zero size,
    {
        if(pF0 != NULL) forcerVecWipe();// discard any existing forcer vectors
        pF0 = new CForcer[fsiz];        // alloc. memory needed for 'f0'
        //... and any other forcer-vectors (see stateVecSizer() for examples)
        if(NULL == pF0)
        {
            cout << "\nERROR: CPartSys::forceVecsizer() failed to alloc memory!?!?! \n" << endl;
            forcerCount = 0;
        }
        else
        {
            forcerCount = fsiz;
        }
    }
    else
    {
        cout << "CPartSys::forcerVecSize(val) val<=0?!? call forcerVecWipe() instead!" << endl;
        forcerVecWipe();
    }
}

void CPartSys::forcerVecWipe(void)
//------------------------------------------------------------------------------
// Discard all currently-allocated memory (if any) for all forcer-vectors.
{
    // Want to expand to more forcer-vectors?
    // Look at CPartSys::stateVecWipe() to see how...
    if(pF0!=NULL)               // Do we have *ANY* CForcer objects
                                // in any forcer-vector?
    {   // Yes.  then we should have the same number of CForcer objects at each
        // and every forcer-vector pointer.  Delete them all, but complain if
        // any are missing:
        if(pF0==NULL)       // make sure vector isn't empty, and
            cout << "CPartSys::forcerVecWipe()--empty f0 vector ?!?" << endl;
        else delete(pF0);   // empty it.
    }
    else
    {
        cout << "\n CPartSys::forceVecWipe() called with all forcer-vectors empty !?! \n" << endl;
    }
}

void CPartSys::forcerVecDraw(CForcer *pF, CPart *pS)
//------------------------------------------------------------------------------
// Use openGL to draw all enabled force-making objects in the forcer-vector
// found at location pF. (e.g. bubbles, springs). You NEED the state-vector
// at pS to find the location of the particles that define the ends of each
// CForcer object that describes a spring.
{
int i,j,tp;

    for(j=0; j<forcerCount; j++)        // for all force-applying objects,
    {
        tp = pF[j].forceType;
        if(tp <= 0) continue;           // ignore disabled force-making objects.
        switch(tp)
        {
            case F_BUBBLE:
                glDisable(GL_LIGHTING); // revert to simple colors;
                glPushMatrix();         // save current matrix;
                    glTranslated(pF[j].bub_ctr[0],
                                 pF[j].bub_ctr[1],
                                 pF[j].bub_ctr[2]); // move to bubble center;
                    // draw a wire sphere, same diameter as bubble;
                    glColor3f(0.3,0.3,0.3);         // make it dim gray.
                    glutWireSphere(pF[j].bub_radius, 8,7);
                glPopMatrix();
                glEnable(GL_LIGHTING);      // restore lighting and shading.
                break;
            case F_SPRING:      // draw a single spring:
                // j-th force-making object's integer e0,e1 values hold the
                // index# of the particles it connects; use this index in the
                // current state s0 to find the particle in the pList array of
                // CPart objects. Within this
                // CPart object, we find position of the particle.
                glDisable(GL_LIGHTING); // revert to simple colors
                glColor3d(0.6, 0.0, 0.0);   // dark red springs
                glBegin(GL_LINES);
                    glVertex3d( pS[ pF[j].e0].val[PART_XPOS ],
                                pS[ pF[j].e0].val[PART_YPOS ],
                                pS[ pF[j].e0].val[PART_ZPOS ]);
                    // other end of spring
                    glVertex3d( pS[ pF[j].e1].val[PART_XPOS ],
                                pS[ pF[j].e1].val[PART_YPOS ],
                                pS[ pF[j].e1].val[PART_ZPOS ]);
                glEnd();
                break;
            case F_SPRINGSET:
                glDisable(GL_LIGHTING); // revert to simple colors
                glColor3d(0.6, 0.0, 0.0);   // dark red springs
                glBegin(GL_LINES);
                for(i=0; i<pF[j].numSprings; i++) // for each spring,
                {
                    glVertex3d(pS[ pF[j].pE0[i] ].val[PART_XPOS],
                               pS[ pF[j].pE0[i] ].val[PART_YPOS],
                               pS[ pF[j].pE0[i] ].val[PART_ZPOS]);
                    glVertex3d(pS[ pF[j].pE1[i] ].val[PART_XPOS],
                               pS[ pF[j].pE1[i] ].val[PART_YPOS],
                               pS[ pF[j].pE1[i] ].val[PART_ZPOS]);
                }
                glEnd();
                break;
            default:
                cout << "CPartSys::forcerVecDraw() -- UNIDENTIFIED force-type ?!?" << endl;
                break;
        }
    }
}

void CPartSys::wallVecSizer(int nuSiz)
//------------------------------------------------------------------------------
// If constraint vectors contain any CWall objects, DELETE them all; make empty
// constraint vectors (e.g. pW0).  Then for each constraint-vector, create a new
// dynamically-alloc'd array that holds 'wallCount' constraints (CWall objects)
// and point to it.
// (Remember, a constraint vector is a pointer to the first object in an array
//   of CWall objects).
{
    if(nuSiz >0)
    {
        if(pC0 != NULL) wallVecWipe();  // discard any existing constraint vectors
        pC0     = new CWall [nuSiz];    // Create new constraint vector
        //... add on more here as needed...
        if(pC0==NULL)                   // if ANY allocation failed,
        {
            cout << "\n!!!CPartSys::wallVecSizer() allocation error!!!\n" << endl;
            exit(1);
        }
        else                            // if it succeeded, record array size;
        {
            wallCount = nuSiz;
        }
    }
    else
    {
        cout << "?!? CPartSys::wallVecSizer(val) with val<=0?!? call wallVecWipe() instead!" << endl;
        wallVecWipe();
    }

}


void CPartSys::wallVecWipe(void)
//------------------------------------------------------------------------------
// Discard all currently-allocated memory (if any) for all wall-vectors
// (e.g. constraint-vectors).
{
    // Want to expand to handle more wall-vectors (constraint-vectors)?
    // Look at CPartSys::stateVecWipe() to see how...
    if(pC0!=NULL)               // Do we have *ANY* CWall objects
                                // in any constraint-vector?
    {   // Yes.  then we should have the same number of CWall objects at each
        // and every constraint-vector pointer.  Delete them all, but complain
        // if any are missing:
        if(pF0==NULL)       // make sure vector isn't empty, and
            cout << "CPartSys::wallVecWipe()--empty w0 vector ?!?" << endl;
        else delete(pC0);   // empty it.
    }
    else
    {
        cout << "\n CPartSys::wallVecWipe() called with all wall-vectors empty !?! \n" << endl;
    }
}

void CPartSys::wallVecDraw(CWall *pW)
//------------------------------------------------------------------------------
// Draw on-screen all CWall constraint-making objects found at the pointer pW.
// (number of CWal objects is CPartSys member 'wallCount'
{
int i;
GLfloat P_diffuse[]   = {0.1, 0.5, 1.0, 1.0};	// light cyan,
GLfloat P_specular[]  = {1.0, 1.0, 1.0, 1.0};  //  white specular,
GLfloat P_shiny = 40.0; // medium shinyness

	glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, P_diffuse);
	glMaterialfv(GL_FRONT, GL_SPECULAR, P_specular);
	glMaterialfv(GL_FRONT, GL_SHININESS, &P_shiny);

    // in whatever openGL coordinate system we're currently using,
    // in whatever glMatrixMode() settings we have (should be GL_MODELVIEW)
    // Draw each constraint;
    for(i=0; i< wallCount; i++)         // for each particle,
    {
        if(pW[i].wallType <= 0) continue;   // skip dead constraints
        glPushMatrix();                     // save current matrix, then
        glEnable(GL_LIGHTING);              // be sure to use our material,

//=========================================
// PUT YOUR CONSTRAINT-DRAWING CODE HERE!
//=========================================

        glPopMatrix();                      // restore the saved matrix.///
    }
}

void CPartSys::applyAllForces(CPart *pS, CForcer *pF)
//------------------------------------------------------------------------------
// Complete the state at pS by finding the total force applied to each particle
// by the list of CForcer objects found at pF.  (pF[i] for 0<= i < forcerCount).
{
int i,j;

    // Before we can find the time-derivative of any given state pS, we must
    // first complete state by finding the total of all forces that impinge on
    // each particle in the state vector pS.
    // Begin by clearing the 'force accumulators' for each particle in ss:

//============================
// WRITE YOUR OWN CODE HERE
//=============================

    // Note that we step through the FORCERs one at a time, not the particles;
    // this requires us to use 'force accumulators' for each particle, BUT it
    // lets define force-making objects that apply to just a few particles
    // such as springs, and selective force-makers that change over time, such
    // as forces for boids that depend on nearest-neighbor particle positions.

    // Accumulate all forces:

//============================
// WRITE YOUR OWN CODE HERE
//=============================

    for(j=0; j< forcerCount; j++)            // for each force-making object,
    {
        switch(pF[j].forceType)    // select what kind of force to apply;
        {
        case -F_GRAV_E: // disabled gravity; do nothing;
            break;
        case F_GRAV_E:      // Earth-gravity: pulls all particles 'downward'.

            for(i=0; i<partCount; i++)    // for every particle,
            {
                if( pS[i].partType <= 0)
                {
                    cout << "\nCPartSys::applyAllForces() F_GRAV_E dead particle:" << i << endl;
                    continue; //skip 'dead' particles
                }             // F=ma: gravity's accel is gravConst, so the
                // magnitude of gravity's force is (mass*gravConst);
//============================
// WRITE YOUR OWN CODE HERE
//=============================
/*                mag = pS[i].val[PART_MASS] * pF[j].gravConst;
                // make a vector: scale our unit vector in the 'down' direction:
                pS[i].val[PART_X_FTOT] += mag * pF[j].downDir[0];
                pS[i].val[PART_Y_FTOT] += mag * pF[j].downDir[1];
                pS[i].val[PART_Z_FTOT] += mag * pF[j].downDir[2];
*/
            }
            break;
        case -F_DRAG: // disabled drag; do nothing;
            break;
        case F_DRAG:        // viscous drag: force = -velocity*K_drag.

//============================
// WRITE YOUR OWN CODE HERE
//=============================

            break;
        case -F_BUBBLE:     // disabled bubble; do nothing.
            break;
        case F_BUBBLE:          // Constant inward force towards centerpoint if
                                // particle is > max_radius away from centerpt.
            for(i=0; i<partCount; i++) // for every particle in the state vector pS,
            {
                if(pS[i].partType <= 0)
                {
                    cout << "\nCPartSys::applyAllForces() F_BUBBLE dead particle: " << i << endl;
                    continue;// skip 'dead' particles
                }                // find vector from bubble center to i-th particle;

//============================
// WRITE YOUR OWN CODE HERE
//=============================

            }
            break;
        case F_SPRING:      // ties together 2 particles; distance sets force

            // e0to1== vector from spring end e0 to spring end e1;
            // (e0,e1 are the indices of the CPart objects connected by spring)
//============================
// WRITE YOUR OWN CODE HERE
//=============================
            // stretched spring: >0; spring force's magnitude:
            // mag = stretch * pF[j].K_spring;
            // Use unit vector(xtmp,ytmp,ztmp) to specify direction; scale it
            // by 'mag' to find the force the spring adds to particle e0:
/*
            pS[pF[j].e0].val[PART_X_FTOT] += mag * e0to1[0];
            pS[pF[j].e0].val[PART_Y_FTOT] += mag * e0to1[1];
            pS[pF[j].e0].val[PART_Z_FTOT] += mag * e0to1[2];
            // equal and opposite force on particle e1.
            pS[pF[j].e1].val[PART_X_FTOT] -= mag * e0to1[0];
            pS[pF[j].e1].val[PART_Y_FTOT] -= mag * e0to1[1];
            pS[pF[j].e1].val[PART_Z_FTOT] -= mag * e0to1[2];
*/
            //---------------------------SPRING DAMPING------------------------
            // Real springs have 'damping' -- internal friction forces that
            // resist spring length-change, proportional to length_velocity.
            // How do we find rate-of-change of length?
            // a)-Find the difference in endpoint velocities: (e0.vel - e1.vel)
            // gives net velocity of e0 measured from e1, as if e1 was fixed.
            // --Find how much net velocity is along the spring's length;
            //  (why discard the velocity vector component perpendicular to the
            //   spring? because those SPIN the spring without length change).
            //      --dot-product of net velocity & e0to1 vector (unit length)
            //          gives velocity magnitude in e0to1 direction;
            //          for e0, mag = e0to1 dot (e0.vel - e1.vel)
            // --Damping force OPPOSES velocity, and is proportional to it:
            //      --for e0 endpoint, damping_force = -mag*K_spring*e0to1;
            // find net velocity............................

//============================
// WRITE YOUR OWN CODE HERE
//=============================
            break;

        case -F_SPRING:     // disabled spring: do nothing.
        case F_SPRINGSET:   // a large set of identical springs; used for
                            // cloth, rubbery objects, water, etc.
        case -F_SPRINGSET:  // disabled spring-set; do nothing.
        case F_GRAV_P:      // planetary-gravity; all masses attract all masses.
        case -F_GRAV_P:     // disabled planetary-gravity: do nothing.
        case F_WIND:        // force on a particle is fcn of 3D position
        case -F_WIND:       // disabled wind: do nothing.
        case F_CHARGE:      // attract/repel by charge and inverse distance
        case -F_CHARGE:     // disabled charge: do nothing
        case F_MOUSE:       // Spring-like connection to the mouse cursor.
        case -F_MOUSE:      // disabled mouse: do nothing.
        case F_NONE:        // Non-existent force: ignore this CForcer object
        default:
            cout << "\n!!!ERROR!!! CPartSys::applyAllForces(): unfinished/undefined force type!?!!?\n" << endl;
            break;
        }
    }
}

void CPartSys::dotMaker(CPart *pDotDest, CPart *pSrc, CForcer *pF)
//------------------------------------------------------------------------------
// Find pDotDest, the time-derivative of state pSrc made by the combined effect
// of all the force-making objects in the array of force-making objects at
// ffIn.
{
int i,j;

    applyAllForces(pSrc, pF);              // Find sum-of- forces currently
                                            // impinging on each particle.
    // We now have a complete state s0; we know net force on every particle.
    // Use it to construct ssDot, our time-derivative of state, by applying
    // common sense and Newton's law.
    // First, set *ALL* the state-variable derivatives s0dot to default value of
    // zero (no change).  (why? so that forgotten variables don't change value)
    for(i=0; i<partCount; i++)
    {
        for(j=0; j< PART_MAXVAR; j++)
        {
            pDotDest[i].val[j] = 0.0;
        }
    }
    // Now set any non-zero values:
/*
    for(i=0; i<partCount; i++)    // for each particle in the state vector ss,
    {
        pDotDest[i].partType = pSrc[i].partType;
                                            // make sure we match alive-ness.
        // Time-derivative of s0 state vector's position? copy s0 velocity;
        pDotDest[i].val[PART_XPOS] = pSrc[i].val[PART_XVEL];
        pDotDest[i].val[PART_YPOS] = pSrc[i].val[PART_YVEL];
        pDotDest[i].val[PART_ZPOS] = pSrc[i].val[PART_ZVEL];
        // Time-derivative of s0 state's velocity? copy s0 acceleration;
        // Newton's Law says F=ma, so a = F/m.  Be SURE we have nonzero mass...
        pDotDest[i].val[PART_XVEL] = pSrc[i].val[PART_X_FTOT] /
                          (NU_EPSILON + pSrc[i].val[PART_MASS]);
        pDotDest[i].val[PART_YVEL] = pSrc[i].val[PART_Y_FTOT] /
                          (NU_EPSILON + pSrc[i].val[PART_MASS]);
        pDotDest[i].val[PART_ZVEL] = pSrc[i].val[PART_Z_FTOT] /
                          (NU_EPSILON + pSrc[i].val[PART_MASS]);
        // time-derivative of s0 state vector 'age' is +1; we increase by 1.0
        // on every timestep:
        pDotDest[i].val[PART_AGE] = 1.0;
        // Lets leave the time-derivative at zero for mass, force-accumulator,
        // charge, and any/all others.

    }
*/
}

void CPartSys::solver(void)
//------------------------------------------------------------------------------
// Use the solver specified by 'solverType' to take the current state s0 and
// its derivative s0dot, and compute the next state s1 (& its derivative s1dot).
// Then swap current state contents (s0, s0dot) with new state (s1, s1dot);
// we're then ready to render current state.
{
    // Find s0dot, the time-derivative of the current state s0, by applying all
    // the force-making objects in our forcer-vector pF0;
    dotMaker(pS0dot,pS0,pF0);

    // Numerical integration:
    // Find the next state s1 by using s0, s0dot, timestep, and if necessary,
    // the dotMaker() fcn in a suitable numerical integration method.  Note
    // we're trying to approximate the true trajectory of s0 through state space
    // but we can only approximate it; the only thing we know with absolute
    // certainty is the current state s0 and the time-derivative s0dot, which
    // we compute from from s0 and the force-making objects.
    // The numerical integration methods we've discussed in class are:
    //  --Euler integration: s1 = s0 + timestep*s0dot
    //  --Midpoint integration, quadratic and arbitrary Runga-Kutta integration,
    // and finally, implicit integration.
    //
    switch(solverType)
    {
        case NU_SOLV_IMPLICIT:          // Implicit integration
        case NU_SOLV_MIDPOINT:          // Midpoint method (see Pixar tutorial)
        case NU_SOLV_ADAMS_BASH:        // Adams-Bashforth explicit integrator
        case NU_SOLV_RUNGEKUTTA:        // Runga-kutta, degree set by solvDegree
        case NU_SOLV_VERLET:            // Verlet semi-implicit integration
        case NU_SOLV_VEL_VERLET:        // 'velocity-Verlet' semi-implicit integ.
            cout << "\n ERROR!  SOLVER NOT DONE YET: Defaulted to Euler Integ.\n" << endl;
        case NU_SOLV_EULER:             // Euler integration
                                        //  (explicit or forward integration)
        default:                        // (also our default solver)
            stateVecAplusBtimesC(pS1, pS0, pS0dot, timeStep);
                                        // s1 = s0 + s0dot*timestep
            break;
    }
    stateVecSwap(&pS0, &pS1);      // now SWAP contents of s0,s1
}

void CPartSys::drawMe(void)
//------------------------------------------------------------------------------
// Render the current state of the particle system;
// Draw the particles in 's0', the forces (if any) in 'f0', and the constraints
// aka 'walls' (if any) in 'w0'.
{
    stateVecDraw(pS0);         // Draw current state.
    forcerVecDraw(pF0,pS0); // now draw all the forcer-vector's CForcer objects
                            // EXCEPT the springs (F_SPRING, F_SPRINGSET).
                            // We need access to BOTH s0 AND CForcer objects to
                            // draw the springs: CForcer tell us which particles
                            // are connected by springs, but only s0 can tell
                            // us the current location of those particles.
    wallVecDraw(pC0);       // Draw current constraints.

}

void CPartSys::initExample1(void)
//------------------------------------------------------------------------------
// All-in-one function to create an interesting example particle system with
// a cube of particles.
{

int nuCount;
int i,j,k, ijkMax, tot;
GLdouble xoff, yoff, zoff;

    ijkMax = 8;
    nuCount = ijkMax * ijkMax * ijkMax;

    if(nuCount < 3) nuCount = 3;        //
    // Set the state vectors, which contain lists of particles.
    stateVecSizer(nuCount);

    // Set up our initial set of particles to form a uniformly-spaced 3D grid
    // of points with one corner at the origing, and with uniform velocity.
    // x0ff,yoff,zoff is the starting point for cube of points
#define JT_DIV 50.0       // sets size of the initial cube of points.
    xoff =0.0; yoff = 0.0; zoff = 0.0;
    xoff -=(GLdouble)(ijkMax-1) * 0.5/(nuCount/JT_DIV); // half-cube offset
    yoff -=(GLdouble)(ijkMax-1) * 0.5/(nuCount/JT_DIV);
    zoff -=(GLdouble)(ijkMax-1) * 0.5/(nuCount/JT_DIV);
    tot = 0;                                // count number of particles done.
    for(k=0; k<ijkMax; k++)                 // for each plane of the cube,
    {
        for(j=0; j<ijkMax; j++)
        {                                   // for each row in that plane,
            for(i=0; (i<ijkMax && tot<nuCount); i++,tot++)
            {                               // for each particle in that row,
                pS0[tot].partType = PTYPE_BALL;       // enable the particle,
                pS0[tot].val[PART_MASS] = 1.0; // set mass,
                // positions : cube-like grid
                pS0[tot].val[PART_XPOS] = xoff + (GLdouble)i* 1.0/(nuCount/JT_DIV);
                pS0[tot].val[PART_YPOS] = yoff + (GLdouble)j* 1.0/(nuCount/JT_DIV);
                pS0[tot].val[PART_ZPOS] = zoff + (GLdouble)k* 1.0/(nuCount/JT_DIV);
                // velocity: all the same
                pS0[tot].val[PART_XVEL] = 0.0;
                pS0[tot].val[PART_YVEL] = 0.0;
                pS0[tot].val[PART_ZVEL] = 0.0;
                // clear the force-accumulator for every particle
                pS0[tot].val[PART_X_FTOT] = 0.0;
                pS0[tot].val[PART_Y_FTOT] = 0.0;
                pS0[tot].val[PART_Z_FTOT] = 0.0;
            }
        }
    }
#undef JT_DIV
    stateVecAssign(pS1, pS0);      // make 'next' state same as 'current' state,
    stateVecZero(pS0dot);          // init all time-derivatives to zero.
    stateVecZero(pS1dot);
    // Set up an interesting set of force-making objects;
    forcerVecSizer(3);    // Make 3 forces: gravity, drag, and a 'bubble' where
j=0;
    // Earth gravity -- pulls in 'downDir' direction
    pF0[j].forceType = F_GRAV_E;        // Earth gravity.
    pF0[j].gravConst  =  3.0;           // accel. due to gravity: 32 ft/sec^2
    pF0[j].downDir[0] =  0.0;           // unit vector in the 'down' direction
    pF0[j].downDir[1] = -1.0;
    pF0[j].downDir[2] =  0.0;           //  gravity pulls in -Z direction.
j=1;
    // Viscous drag: force = -velocity*K_drag
    pF0[j].forceType = F_DRAG;
    pF0[j].K_drag = 0.6;                // damping strength.
j=2;
    // 'Bubble' within an inward force field:
    pF0[j].forceType = F_BUBBLE;
    pF0[j].bub_radius = 1.0; // 1.0;    // anything outside this radius,
    pF0[j].bub_ctr[0] = 0.0;            // measured from this center,
    pF0[j].bub_ctr[1] = 0.0;
    pF0[j].bub_ctr[2] = 0.0;
    pF0[j].bub_force = -5.0;            // gets this force towards bub_ctr
    runNotStep = false;                 // show initial frame only; user must
                                        // hit 'P' key to run, spacebar to  step

}

void CPartSys::initExample2(void)
//------------------------------------------------------------------------------
// All-in-one function to create an interesting example particle system.
{

//============================
// WRITE YOUR OWN CODE HERE
//=============================

}

void CPartSys::initExample3(void)
//------------------------------------------------------------------------------
// All-in-one function to create an interesting example particle system.
// Makes a simple tetrahedron with springs, in a bubble, with gravity & drags.
{

int nuCount;
int i,j;
GLdouble xoff, yoff, zoff;

    nuCount = 4;                      // make a spring-mass tetrahedron
    stateVecSizer(nuCount);

    xoff = -0.3;    yoff = -0.3;    zoff = -0.3;
// use some macros to avoid tedious re-typing of almost-identical statements:
#define NU_X(v)  (pS0[i].val[PART_XPOS] = xoff + v)
#define NU_Y(v)  (pS0[i].val[PART_YPOS] = yoff + v)
#define NU_Z(v)  (pS0[i].val[PART_ZPOS] = zoff + v)
    i=0; NU_X(0.0); NU_Y(0.0); NU_Z(0.0);       // particle 0 at origin,
    i=1; NU_X(0.6); NU_Y(0.0); NU_Z(0.0);       //          1 at +x
    i=2; NU_X(0.0); NU_Y(0.6); NU_Z(0.0);       //          2 at +y
    i=3; NU_X(0.0); NU_Y(0.0); NU_Z(0.6);       //          3 at +z.
#undef NU_X   // (locals only!)
#undef NU_Y
#undef NU_Z
    for(i=0; i<nuCount; i++)    // otherwise, the particles are identical:
    {                               // for each particle,
        pS0[i].partType = PTYPE_BALL;     // enable the particle,
        pS0[i].val[PART_MASS] = 0.5;        // set mass,
        // velocity: all the same
        pS0[i].val[PART_XVEL] = 0.0;
        pS0[i].val[PART_YVEL] = 0.0;
        pS0[i].val[PART_ZVEL] = 0.0;
        // clear the force-accumulator:
        pS0[i].val[PART_X_FTOT] = 0.0;
        pS0[i].val[PART_Y_FTOT] = 0.0;
        pS0[i].val[PART_Z_FTOT] = 0.0;
    }

    stateVecAssign(pS1, pS0);   // make 'next' state same as 'current' state,
    stateVecZero(pS0dot);       // init all time-derivatives to zero.
    stateVecZero(pS1dot);
    // Now set up force-making objects:
    forcerVecSizer(9);          // Make 9 forces: 6 springs for the tetrahedron,
                                // plus gravity, drag, and 1 'bubble' object:
    j=0;
    // Earth gravity -- pulls in 'downDir' direction
    pF0[j].forceType = F_GRAV_E;        // Earth gravity.
    pF0[j].gravConst  =  3.0;           // accel. due to gravity: 32 ft/sec^2
    pF0[j].downDir[0] =  0.0;           // unit vector in the 'down' direction
    pF0[j].downDir[1] = -1.0;
    pF0[j].downDir[2] =  0.0;           //  gravity pulls in -Z direction.
    // Viscous drag: force = -velocity*K_drag
    j=1;
    pF0[j].forceType = F_DRAG;
    pF0[j].K_drag = 0.6;                // damping strength.
    // 'Bubble' within an inward force field:
    j=2;
    pF0[j].forceType = F_BUBBLE;
    pF0[j].bub_radius = 1.0;            // anything outside this radius,
    pF0[j].bub_ctr[0] = 0.0;            // measured from this center,
    pF0[j].bub_ctr[1] = 0.0;
    pF0[j].bub_ctr[2] = 0.0;
    pF0[j].bub_force = -5.0;            // gets this force towards bub_ctr

    //--------------------Springs

    for(j=3; j<9; j++)  // last 6 force-making objects are springs for our
    {                   // tetrahedron.  All these springs have:
        pF0[j].forceType = F_SPRING;
        pF0[j].K_spring = 50.0;
        pF0[j].K_restLength = 0.4;
        pF0[j].K_springdamp = 0.8;
    }       // Now connect each spring to its corresponding particles;

    j=3;    // spring from particle 0 to 1
    pF0[j].e0 = 0;   pF0[j].e1 = 1;
    j=4;    // spring from particle 0 to 2
    pF0[j].e0 = 0;   pF0[j].e1 = 2;
    j=5;    // spring from particle 0 to 3
    pF0[j].e0 = 0;   pF0[j].e1 = 3;
    j=6;    // spring from particle 1 to 2
    pF0[j].e0 = 1;   pF0[j].e1 = 2;
    j=7;    // spring from particle 2 to 3
    pF0[j].e0 = 2;   pF0[j].e1 = 3;
    j=8;    // spring from particle 3 to 1
    pF0[j].e0 = 3;   pF0[j].e1 = 1;

    runNotStep = false;                 // show initial frame only; user must
                                       // hit 'P' key to run, spacebar to  step
}
